Równania liniowej sprężystości
W module tym przedstawimy sformułowanie silne i słabe dla równań liniowej teorii sprężystości, bardzo powszechnie używanych podczas obliczeń za pomocą metody elementów skończonych naprężeń w różnych konstrukcjach.
Rozważmy dla przykładu trójwymiarową platformę podpartą pięcioma podporami stojącą na płaszczyźnie OXY, przedstawioną na rysunku 1. Platforma ta posiada swój ciężar i pojawiają się na niej pewne naprężenia wynikające z oddziaływania pola grawitacyjnego ziemi.
Dla wygodniejszego sformułowania równań różniczkowych cząstkowych opisujących przemieszczenie platformy, będziemy teraz używać oznaczenia \( (x_1,x_2,x_3) \) w dowolnym punkcie platformy.
Poszukiwanym rozwiązaniem będzie tym razem pole wektorowe przemieszczeń
\( \Omega \ni (x_1,x_2,x_3) \rightarrow (u_1(x_1,x_2,x_3),u_2(x_1,x_2,x_3),u_3(x_1,x_2,x_3)) \)
Oznacza to w szczególnosci że przemieszczenia punktów platformy pod jej ciężarem, w kierunku osi \( z \), wynoszą \( u_3(x_1,x_2,x_3) \) w punkcie \( (x_1,x_2,x_3) \). Oczywiście pole grawitacyjne działa jedynie w kierunku osi Z, ale punkty platformy są ze sobą "sklejone" i jeśli punkt platformy przemieszcza się w kierunku osi Z pod wpływem ciężaru platformy, często związane jest to z przesunięciem punktu w kierunkach Y oraz Z, co opisują składowe \( u_1(x_1,x_2,x_3) \) oraz \( u_2(x_1,x_2,x_3) \). Oczywiście przewidzenie jak będą odkształcać się te punkty nie jest rzeczą prostą, i w szczególności wymaga rozwiązania skomplikowanych równań liniowej teorii sprężystości.
Dla uproszczenia zapisu oraz ze względu na stosowanie tradycyjnego zapisu wprowadzamy oznaczenie na pochodne cząstkowe poszczególnych współrzędnych wektora przemieszczeń
\( u_{i,j}=\frac{\partial u_i}{\partial x_j } \)
oznacza to pochodną w kierunku \( j \)-tym ze składowej pola wektorowego \( i \)-tego.
W celu wyprowadzenia równań liniowej teorii sprężystości wprowadza się tak zwany tensor prędkości odkształceń
\( \epsilon_{ij}=\frac{\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}}{2}=\frac{u_{i,j}+u_{j,i } }{2 } \)
oraz tensor naprężeń
\( \sigma_{ij}=c_{ijkl}\epsilon_{kl} \)
Używamy tutaj tak zwanej konwencji sumacyjnej Einsteina, w której zduplikowanie indeksów we wzorze oznacza iż sumujemy po tych indeksach czyli na przykład w równaniu spreżystości
\( \sigma_{ij}=c_{ijkl}\epsilon_{kl } \)
powtórzenie indeksów \( k,l \) oznacza sumowanie po tych indeksach
\( \sigma_{ij}=\sum_{k=1,...,3;l=1,...,3} c_{ijkl}\epsilon_{kl } \).
Mamy tutaj czterowymiarową macierz współczynników \( c_{ijkl},i=1,...,3;j=1,...,3,k=1,...,3;l=1,..,3 \). Określenie tych 3*3*3*3 liczb definiuje mi jednoznacznie materiał z którego zbudowana jest analizowana przeze mnie konstrukcja (w naszym przykładzie platforma). Większość tych współczynników na szczęście będzie równa zero.
Tak zwane sformułowanie silne problemu liniowej teorii sprężystości formułuje się następująco:
\( \sigma_{ij,j}+f_i=0 \) dla i=1,2,3, w obszarze \( \Omega \). Wyraz \( f_i \) oznacza 3 współrzędne wektora siły działającej na naszą konstrukcję, na przykłąd będzie to siła grawitacji \( f=(g_1(x_1,x_2,x_3),g_2(x_1,x_2,x_3),g_3(x_1,x_2,x_3)) \)
Wprowadzamy tutaj następującą konwencję. Indeksy powtórzone po przecinku \( ij,j \) oznaczają obliczanie pochodnej pochodne cząstkowej po indeksowanej zmiennej.
Zapis
\( \sigma_{ij,j}=\frac{\partial \sigma_{ij} }{\partial x_j } \) oznacza więc sumę pochodnych cząstkowych z tensora naprężeń. Ponadto mamy defacto trzy równania (co ma sens ponieważ szukamy trzech współrzędnych wektora przemieszczeń).
\( \sum_{j=1,2,3}\frac{\partial \sigma_{1j}}{\partial x_j}+f_1=0 \\ \sum_{j=1,2,3}\frac{\partial \sigma_{2j}}{\partial x_j }+f_2=0 \\ \sum_{j=1,2,3}\frac{\partial \sigma_{3j}}{\partial x_j }+f_3=0 \)
czyli uwzględniając defincje tensora naprężeń
\( \sum_{j=1,2,3}\frac{\partial (\sum_{k=1,...,3;l=1,...,3} c_{1jkl}\epsilon_{kl})}{\partial x_j}+f_1=0 \\ \sum_{j=1,2,3}\frac{\partial (\sum_{k=1,...,3;l=1,...,3} c_{2jkl}\epsilon_{kl})}{\partial x_j}+f_2=0 \\ \sum_{j=1,2,3}\frac{\partial (\sum_{k=1,...,3;l=1,...,3} c_{3jkl}\epsilon_{kl})}{\partial x_j }+f_3=0 \)
czyli uwzględniając definicje tensora prędkości odkształceń
\( \sum_{j=1,2,3}\frac{\partial (\sum_{k=1,...,3;l=1,...,3} c_{1jkl}(\frac{u_{k,l}+u_{l,k}}{2}))}{\partial x_j }+f_1=0 \\ \sum_{j=1,2,3}\frac{\partial (\sum_{k=1,...,3;l=1,...,3} c_{2jkl}(\frac{u_{k,l}+u_{l,k}}{2}))}{\partial x_j }+f_2=0 \\ \sum_{j=1,2,3}\frac{\partial (\sum_{k=1,...,3;l=1,...,3} c_{3jkl}\frac{u_{k,l}+u_{l,k}}{2})}{\partial x_j }+f_3=0 \)
Nasze niewiadome tutaj to pole wektorowe przemieszczeń: Szukamy
\( \Omega \ni (x_1,x_2,x_3) \rightarrow (u_1(x_1,x_2,x_3),u_2(x_1,x_2,x_3),u_3(x_1,x_2,x_3)) \)
Żeby rozwiązać równania liniowej teorii sprężystości na danej konstrukcji, konieczne jest zdefiniowanie warunków brzegowych, czyli określenie tego co się dzieje na brzegu konstrukcji. Na części brzegu możemy wprowadzić warunek Dirichleta, na przykład określając że podstawa konstrukcji jest przymocowana do podłoża, i dlatego przemieszczenie podstawy konstrukcji jest zerowe
\( u_1(x_1,x_2,x_3)=0, u_2(x_1,x_2,x_3)=0, u_3(x_1,x_2,x_3)=0, \textrm{ dla } (x_1,x_2,x_3)\in \Gamma_D \)
Możliwe jest również określenie warunku brzegowego Neumanna
\( \sigma_{ij}n_j = h_i \textrm{ dla } (x_1,x_2,x_3)\in \Gamma_N \)
gdzie \( n(x_1,x_2,x_3)=(n_1(x_1,x_2,x_3),n_2(x_1,x_2,x_3),n_3(x_1,x_2,x_3)) \) oznacza wersor prostopadły do brzegu w punkcie \( (x_1,x_2,x_3) \), oraz \( h=(h_1(x_1,x_2,x_3),h_2(x_1,x_2,x_3),h_3(x_1,x_2,x_3)) \) oznacza zadaną funkcję na brzegu Neumanna \( \Gamma_N \), czyli uwzględniając konwencje sumacyjną Einsteina
\( \sigma_{11}n_2+\sigma_{12}n_2+\sigma_{13}n_3 = h_1 \\ \sigma_{21}n_2+\sigma_{22}n_2+\sigma_{23}n_3 = h_2 \\ \sigma_{31}n_2+\sigma_{32}n_2+\sigma_{33}n_3 = h_3 \)
czyli uwzględniając definicje tensora naprężeń
\( (\sum_{k=1,...,3;l=1,...,3} c_{11kl}\epsilon_{kl})n_2+ (\sum_{k=1,...,3;l=1,...,3} c_{12kl}\epsilon_{kl})n_2+ (\sum_{k=1,...,3;l=1,...,3} c_{13kl}\epsilon_{kl})n_3 = h_1 \\ (\sum_{k=1,...,3;l=1,...,3} c_{21kl}\epsilon_{kl})n_2+(\sum_{k=1,...,3;l=1,...,3} c_{22kl}\epsilon_{kl})n_2+ (\sum_{k=1,...,3;l=1,...,3} c_{23kl}\epsilon_{kl})n_3 = h_2 \\ (\sum_{k=1,...,3;l=1,...,3} c_{31kl}\epsilon_{kl})n_2+ (\sum_{k=1,...,3;l=1,...,3} c_{32kl}\epsilon_{kl})n_2+ (\sum_{k=1,...,3;l=1,...,3} c_{33kl}\epsilon_{kl })n_3 = h_3 \).
Warunek Neumanna określa siły w równaniach teorii sprężystości
Sformułowanie słabe (zwane też wariacyjnym) które nadawać się będzie do implementacji za pomocą metody elementów skończonych podajemy na podstawie książki
[1]
Sformułowanie słabe (wariacyjne) podane poniżej obowiązuje dla pewnej szerokiej klasy materiałów dla której własności mechaniczne da się opisać za pomocą dwóch parametrów, modułu Younga \( E \) oraz współczynnika Poissona \( \mu \).
Dla zadanego pola sił
\( \Omega \ni (x_1,x_2,x_3) \rightarrow f=(g_1(x_1,x_2,x_3),g_2(x_1,x_2,x_3),g_3(x_1,x_2,x_3)) \)
oraz zadanej funkcji warunku brzegowego Neumanna
\( \Omega \ni (x_1,x_2,x_3) \rightarrow h=(h_1(x_1,x_2,x_3),h_2(x_1,x_2,x_3),h_3(x_1,x_2,x_3)) \)
oraz zadanego tensora (macierzy czterowymiarowej) współczynników materiałowych \( c_{ijkl} \) obliczyć pole przemieszczeń materiału konstrukcji
\( \Omega \ni (x_1,x_2,x_3) \rightarrow (u_1(x_1,x_2,x_3),u_2(x_1,x_2,x_3),u_3(x_1,x_2,x_3)) \)
spełniające równania liniowej sprężystości
\( a(w,u)=(w,f)+(w,h)_{\Gamma} \)
dla dowolnych funkcji testujących \( w=(w_1(x_1,x_2,x_3),w_2(x_1,x_2,x_3),w_3(x_1,x_2,x_3)) \),
gdzie
\( (w,f)=\int_{\Omega}w_1(x_1,x_2,x_3)f_1(x_1,x_2,x_3)dx+\int_{\Omega}w_2(x_1,x_2,x_3)f_2(x_1,x_2,x_3)dy+\int_{\Omega}w_3(x_1,x_2,x_3)f_3(x_1,x_2,x_3)dz \\ (w,h)=\int_{\Gamma}w_1(x_1,x_2,x_3)h_1(x_1,x_2,x_3)dS+\int_{\Gamma}w_2(x_1,x_2,x_3)h_2(x_1,x_2,x_3)dS+\int_{\Gamma}w_3(x_1,x_2,x_3)f_3(x_1,x_2,x_3)dS \\a(w,u)=\int_{\Omega}\epsilon(w)^T D \epsilon(u) dxdydz \),
gdzie
\( \epsilon(u) = \begin{bmatrix} u_{1,1} \\ u_{2,2} \\ u_{3,3} \\ u_{2,3}+u_{3,2} \\ u_{1,3}+u_{3,1} \\ u_{1,2}+u_{2,1} \\ \end{bmatrix}, \epsilon(w) = \begin{bmatrix} w_{1,1} \\ w_{2,2} \\ w_{3,3} \\ w_{2,3}+w_{3,2} \\ w_{1,3}+w_{3,1} \\ w_{1,2}+w_{2,1} \\ \end{bmatrix} \)
natomiast
\( D = \frac{E}{(1+2\mu)(1-2\mu) } \begin{bmatrix} 1-\mu & \mu & \mu & 0 & 0 & 0 \\ \mu & 1-\mu & \mu & 0 & 0 & 0 \\ \mu & \mu & 1-\mu & 0 & 0 & 0 \\ 0 & 0 & 0& \frac{1-2\mu}{2} & 0 & 0 \\ 0 & 0 & 0& 0 & \frac{1-2\mu}{2} & 0 \\ 0 & 0 & 0& 0 & 0& \frac{1-2\mu}{2} \\ \end{bmatrix} \)
wyprowadzona jest przy założeniu że ogólny tensor własności materiału \( c_{ijkl} \) da się opisać za pomocą dwóch parametrów modułu Younga \( E \) oraz współczynnika Poissona \( \mu \). Podane sformułowanie wariacyjne jest stosunkowo wygodne do implementacji. Można w tym celu zastosować izogeometryczną metodę elementów skończonych, lub klasyczną metodę elementów skończonych.
Rozważmy następujący przykład.
Przyjmijmy że nasz obszar \( \Omega \) to platforma o wymiarach 16 na 16 bloków, podparta pięcioma podporamiumieszczonymi pod blokami (1,1), (16,1), (1,16), (16,16) oraz (8,8), tak jak pokazano to na rysunku Rys. 1.
Przyjmujemy zerowy warunek Dirichleta u podstawy platformy (czyli zakładamy że nóżki są przyczepione do płaskiego stabilnego podłoża), natomiast na reszcie brzegu nie przyjmujemy żadnego warunku brzegowego. Takie założenie jest wystarczające ponieważ umocowanie platformy u podstaw nóżek daje nam możliwosć uzyskania jednoznacznego rozwiązania (obliczenia odkształceń wynikających z ciężaru platofrmy). Zakładamy siłe grawitacyjną \( f=(0,0,-1) \), oraz \( h=0 \) (brak warunku brzegowego Neumanna.
Mamy więc
\( a(w,u)=(w,f)+(w,h)_{\Gamma} \\ a(w,u)=\int_{\Omega}\epsilon(w)^T D \epsilon(u) dxdydz \\(w,f)=\int_{\Omega}w_1(x_1,x_2,x_3)f_1(x_1,x_2,x_3)dx+\int_{\Omega}w_2(x_1,x_2,x_3)f_2(x_1,x_2,x_3)dy+\int_{\Omega}w_3(x_1,x_2,x_3)f_3(x_1,x_2,x_3)dz= \\ \int_{\Omega}0*f_1(x_1,x_2,x_3)dx+\int_{\Omega}0*f_2(x_1,x_2,x_3)dy+\int_{\Omega}(-1)*f_3(x_1,x_2,x_3)dz \\ (w,h)=(w,0)=0 \)
Ponadto przyjmujemy moduł Younga \( E=1[GPa] \)
oraz \( \mu=0.3 \) czyli
\( D = \frac{1}{(1.6)(0.4) } \begin{bmatrix} 0.7 & 0.3 & 0.3 & 0 & 0 & 0 \\ 0.3 & 0.7 & 0.3 & 0 & 0 & 0 \\ 0.3 & 0.3 & 0.7 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.2 & 0 & 0\\ 0 & 0 & 0& 0 & 0.2 & 0 \\ 0 & 0 & 0 & 0 & 0& 0.2 \\ \end{bmatrix} [GPa] \)
Wówczas powierzchnię platformy oraz każdą z pięciu podstaw potraktować można jako osobną grupę elementów i opisać za pomocą następujących wektorów węzłów
Platforma:
[0 0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 16 16] x [0 0 0 1 1 1] x [0 0 0 1 1 1]
Podstawy:
[0 0 0 1 1 1] x [0 0 0 1 1 1] x [0 0 0 10 10 10]
[15 15 15 16 16 16] x [15 15 15 16 16 16] x [0 0 0 10 10 10]
[0 0 0 1 1 1] x [15 15 15 16 16 16] x [0 0 0 10 10 10]
[15 15 15 16 16 16] x [0 0 0 1 1 1] x [0 0 0 10 10 10]
[7 7 7 8 8 8] x [7 7 7 8 8 8] x [0 0 0 10 10 10]
Każdy z komponentów wektorowego pola odkształceń przybliżamy za pomocą kombinacji liniowej funkcji B-spline na poszczególnych grupach elementów
\( u^1(x_1,x_2,x_3) = \sum_{i,j,k=1,...,N_x,N_y,N_z} u^1_{i,j,k} B^{x_1}_{i,p}(x_1)B^{x_2}_{j,p}(x_2)B^{x_3}_{k,p }(x_3) \\ u^2(x_1,x_2,x_3) = \sum_{i,j,k=1,...,N_x,N_y,N_z} u^2_{i,j,k} B^{x_1}_{i,p}(x_1)B^{x_2}_{j,p}(x_2)B^{x_3}_{k,p }(x_3) \\ u^3(x_1,x_2,x_3) = \sum_{i,j,k=1,...,N_x,N_y,N_z} u^3_{i,j,k} B^{x_1}_{i,p}(x_1)B^{x_2}_{j,p}(x_2)B^{x_3}_{k,p }(x_3) \)
W naszym problemie liniowej sprężystości rozwiązujemy równocześnie trzy składowe pola przemieszczeń i wkomponowują się one razem w jedną globalną macierz.
Wstawiamy dyskretyzacje do sformułowania słabego, oraz za funkcje testujące \( w \) również przyjmujemy funkcje B-spline, generujemy układ równań, rozwiązujemy go solwerem dokładnym i uzyskujemy rozwiązanie przedstawione na rysunku Rys. 2. Rysunek przedstawią normę z odkształceń, zdefiniowaną jako \( \| u^1(x_1,x_2,x_3)^2+u^2(x_1,x_2,x_3)^2+u^3(x_1,x_2,x_3)^2 \| \). Uzyskane rozwiązanie nie jest symetryczne, ponieważ platforma posiada oprócz czterech podpór umieszczonych w rogach, również piątą podporę umieszczoną niesymetrycznie.